home *** CD-ROM | disk | FTP | other *** search
/ Amiga Collections: Camelot / Camelot 098 (1990-12)(Swedish User Group of Amiga)(SE)(PD)[WB].zip / Camelot 098 (1990-12)(Swedish User Group of Amiga)(SE)(PD)[WB].adf / XLisp-Stat / statistics.lsp < prev    next >
Lisp/Scheme  |  1990-10-03  |  11KB  |  344 lines
  1. ;;;;
  2. ;;;; statistics.lsp XLISP-STAT statistics functions
  3. ;;;; XLISP-STAT 2.1 Copyright (c) 1988, by Luke Tierney
  4. ;;;; Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz
  5. ;;;; You may give out copies of this software; for conditions see the file
  6. ;;;; COPYING included with this distribution.
  7. ;;;;
  8.  
  9. (provide "statistics")
  10.  
  11. ;;;;
  12. ;;;; Data File Reading 
  13. ;;;;
  14.  
  15. (defun count-file-columns (fname)
  16. "Args: (fname)
  17. Returns the number of lisp items on the first nonblank line of file FNAME."
  18.   (with-open-file (f fname)
  19.     (if f
  20.         (let ((line (do ((line (read-line f) (read-line f))) 
  21.                         ((or (null line) (< 0 (length line))) line))))
  22.           (if line
  23.               (with-input-from-string (s line)
  24.                    (do ((n 0 (+ n 1)) (eof (gensym))) 
  25.                     ((eq eof (read s nil eof)) n))))))))
  26.  
  27. (defvar *xlisptable* *readtable*)
  28.  
  29. #+(and (not macintosh) (not amiga) dialogs)
  30. (defun open-file-dialog (&optional set) 
  31.   (get-string-dialog "Enter a data file name:"))
  32. #-dialogs
  33. (defun open-file-dialog (&optional set) 
  34.   (error "You must provide a file name explicitly"))
  35.  
  36. (defun read-data-file (&optional (file (open-file-dialog t)))
  37. "Args:  (file)
  38. Returns a list of all lisp objects in FILE. FILE can be a string or a symbol,
  39. in which case the symbol'f print name is used."
  40.   (if file
  41.       (let ((oldtable *readtable*)
  42.             (oldbreak *breakenable*)
  43.             (eof (gensym)))
  44.         (setq *readtable* *xlisptable*)
  45.         (setq *breakenable* nil)
  46.         (with-open-file (f file)
  47.           (if f
  48.               (unwind-protect
  49.                (do* ((r (read f nil eof) (read f nil eof))
  50.                      (x (list nil))
  51.                      (tail x (cdr tail)))
  52.                     ((eq r eof) (cdr x))
  53.                     (setf (cdr tail) (list r)))
  54.                (setq *breakenable* oldbreak)
  55.                (setq *readtable* oldtable)))))))
  56.  
  57. ;;; New definition to avoid stack size limit in apply
  58. (defun read-data-columns (&optional (file (open-file-dialog t))
  59.                                     (cols (if file 
  60.                                               (count-file-columns file))))
  61. "Args: (&optional file cols)
  62. Reads the data in FILE as COLS columns and returns a list of lists representing the columns."
  63.   (if (and file cols)
  64.       (transpose (split-list (read-data-file file) cols))))
  65.  
  66. #+(or unix amiga)
  67. (defun load-data (file)
  68. "Args: (file)
  69. Read in data file from the data examples library."
  70.   (if (load (format nil "~aData/~a" *default-path* file))
  71.       t
  72.       (load (format nil "~aExamples/~a" *default-path* file))))
  73.  
  74. #+(or unix amiga)
  75. (defun load-example (file)
  76. "Args: (file)
  77. Read in lisp example file from the examples library."
  78.   (if (load (format nil "~aExamples/~a" *default-path* file))
  79.       t
  80.       (load (format nil "~aData/~a" *default-path* file))))
  81. #+macintosh
  82. (defun load-data (s) (require s (concatenate 'string ":Data:" s)))
  83. #+macintosh
  84. (defun load-example (s) (require s (concatenate 'string ":Examples:" s)))
  85.  
  86. ;;;;
  87. ;;;; Listing and Saving Variables and Functions
  88. ;;;;
  89.  
  90. (defvar *variables*)
  91. (defvar *ask-on-redefine*)
  92.  
  93. (defmacro def (symbol value)
  94. "Syntax: (def var form)
  95. VAR is not evaluated and must be a symbol.  Assigns the value of FORM to
  96. VAR and adds VAR to the list *VARIABLES* of def'ed variables. Returns VAR.
  97. If VAR is already bound and the global variable *ASK-ON-REDEFINE*
  98. is not nil then you are asked if you want to redefine the variable."
  99.   `(unless (and *ask-on-redefine*
  100.                 (boundp ',symbol)
  101.                 (not (y-or-n-p "Variable has a value. Redefine?")))
  102.            (pushnew ',symbol *variables*)
  103.            (setf ,symbol ,value)
  104.            ',symbol))
  105.   
  106. (defun variables-list () 
  107.     (mapcar #'intern (sort-data (mapcar #'string *variables*))))
  108.  
  109. (defun variables ()
  110. "Args:()
  111. Returns a list of the names of all def'ed variables to STREAM"
  112.   (if *variables*
  113.       (mapcar #'intern (sort-data (mapcar #'string *variables*)))))
  114.   
  115. (defun savevar (vars file)
  116. "Args: (vars file-name-root)
  117. VARS is a symbol or a list of symbols. FILE-NAME-ROOT is a string (or a symbol
  118. whose print name is used) not endinf in .lsp. The VARS and their current values
  119. are written to the file FILE-NAME-ROOT.lsp in a form suitable for use with the
  120. load command."
  121.   (let ((f (open (strcat (string file) ".lsp") :direction :output))
  122.         (vars (if (consp vars) vars (list vars)))
  123.         (oldbreak *breakenable*))
  124.     (setq *breakenable* nil)
  125.     (unwind-protect
  126.       (mapcar
  127.         (lambda (x)
  128.             (if (objectp (symbol-value x))
  129.                 (print `(def ,x ,(send (symbol-value x) :save)) f)
  130.                 (print `(def ,x ',(symbol-value x)) f)))
  131.         vars)
  132.       (setq *breakenable* oldbreak)
  133.       (close f))
  134.     vars))
  135.  
  136. (defun undef (v)
  137. "Args: (v)
  138. If V is the symbol of a defined variable the variable it is unbound and
  139. removed from the list of defined variables. If V is a list of variable
  140. names each is unbound and removed. Returns V."
  141.   (dolist (s (if (listp v) v (list v)))
  142.           (when (member s *variables*)
  143.                 (setq *variables* (delete s *variables*))
  144.                 (makunbound s)))
  145.   v)
  146.         
  147.  
  148. ;;;;
  149. ;;;; Basic Summary Statistics
  150. ;;;;
  151.  
  152. (defun standard-deviation (x)
  153. "Args: (x)
  154. Returns the standard deviation of the elements x. Vector reducing."
  155.   (let ((n (count-elements x))
  156.         (r (- x (mean x))))
  157.     (sqrt (* (mean (* r r)) (/ n (- n 1))))))
  158.  
  159. (defun quantile (x p)
  160. "Args: (x p)
  161. Returns the P-th quantile(s) of sequence X. P can be a number or a sequence."
  162.   (let* ((x (sort-data x))
  163.          (n (length x))
  164.          (np (* p (- n 1)))
  165.          (low (floor np))
  166.          (high (ceiling np)))
  167.     (/ (+ (select x low) (select x high)) 2)))
  168.     
  169. (defun median (x) 
  170. "Args: (x)
  171. Returns the median of the elements of X."
  172.   (quantile x 0.5))
  173.  
  174. (defun interquartile-range (x) 
  175. "Args: (number-data)
  176. Returns the interquartile range of the elements of X."
  177.   (apply #'- (quantile x '(0.75 0.25))))
  178.  
  179. (defun fivnum (x) 
  180. "Args: (number-data)
  181. Returns the five number summary (min, 1st quartile, medinan, 3rd quartile,
  182. max) of the elements X."
  183.   (quantile x '(0 .25 .5 .75 1)))
  184.  
  185. (defun covariance-matrix (&rest args)
  186. "Args: (&rest args)
  187. Returns the sample covariance matrix of the data columns in ARGS. ARGS may
  188. consist of lists, vectors or matrices."
  189.   (let ((columns (apply #'append 
  190.                         (mapcar (lambda (x) 
  191.                                   (if (matrixp x) (column-list x) (list x)))
  192.                                 args))))
  193.     (/ (cross-product (apply #'bind-columns 
  194.                              (- columns (mapcar #'mean columns))))
  195.        (- (length (car columns)) 1))))
  196.  
  197. ;;;;
  198. ;;;; Basic Sequence Operations
  199. ;;;;
  200.  
  201. (defun difference (x)
  202. "Args: (x)
  203. Returns differences for a sequence X."
  204.   (let ((n (length x)))
  205.     (- (select x (iseq 1 (1- n))) (select x (iseq 0 (- n 2))))))
  206.  
  207. (defun rseq (a b num)
  208. "Args: (a b num)
  209. Returns a list of NUM equally spaced points starting at A and ending at B."
  210.   (+ a (* (iseq 0 (1- num)) (/ (- b a) (1- num)))))
  211.  
  212.  
  213. ;;;;
  214. ;;;; Linear Algebra Functions
  215. ;;;;
  216.  
  217. (defun matrix (dim data)
  218. "Args: (dim data)
  219. returns a matrix of dimensions DIM initialized using sequence DATA
  220. in row major order." 
  221.   (let ((dim (coerce dim 'list))
  222.         (data (coerce data 'list)))
  223.     (make-array dim :initial-contents (split-list data (nth 1 dim)))))
  224.  
  225. (defun print-matrix (a &optional (stream *standard-output*))
  226. "Args: (matrix &optional stream)
  227. Prints MATRIX to STREAM in a nice form that is still machine readable"
  228.   (unless (matrixp a) (error "not a matrix - ~a" a))
  229.   (let ((size (min 15 (max (map-elements #'flatsize a)))))
  230.     (format stream "#2a(~%")
  231.     (dolist (x (row-list a))
  232.             (format stream "    (")
  233.             (let ((n (length x)))
  234.               (dotimes (i n)
  235.                        (let ((y (aref x i)))
  236.                          (cond
  237.                            ((integerp y) (format stream "~vd" size y))
  238.                            ((floatp y) (format stream "~vg" size y))
  239.                            (t (format stream "~va" size y))))
  240.                        (if (< i (- n 1)) (format stream " "))))
  241.             (format stream ")~%"))
  242.     (format stream "   )~%")
  243.     nil))
  244.  
  245. (defun solve (a b)
  246. "Args: (a b)
  247. Solves A x = B using LU decomposition and backsolving. B can be a sequence
  248. or a matrix."
  249.   (let ((lu (lu-decomp a)))
  250.     (if (matrixp b)
  251.         (apply #'bind-columns
  252.                (mapcar #'(lambda (x) (lu-solve lu x)) (column-list b)))
  253.         (lu-solve lu b))))
  254.         
  255. (defun backsolve (a b)
  256. "Args: (a b)
  257. Solves A x = B by backsolving, assuming A is upper triangular. B must be a
  258. sequence. For use with qr-decomp."
  259.   (let* ((n (length b))
  260.          (sol (make-array n)))
  261.     (dotimes (i n)
  262.              (let* ((k (- n i 1))
  263.                     (val (elt b k)))
  264.                (dotimes (j i)
  265.                         (let ((l (- n j 1)))
  266.                           (setq val (- val (* (aref sol l) (aref a k l))))))
  267.                (setf (aref sol k) (/ val (aref a k k)))))
  268.     (if (listp b) (coerce sol 'list) sol)))
  269. #|
  270. (defun eigenvalues (a)
  271. "Args: (a)
  272. Returns the singular values of A (the eigen values if A is square and
  273. symmetric"
  274.   (coerce (nth 1 (sv-decomp a)) 'list))
  275.  
  276. (defun eigenvectors (a)
  277. "Args: (a)
  278. Returns the left singular vectors of A (the eigen vectors if A is square and
  279. symmetric"
  280.   (column-list (nth 0 (sv-decomp a))))
  281.  
  282. (defun eigen (a)
  283. "Args: (a)
  284. Returns a list of the list of singular values and a list of left singular
  285. vectors of A (the eigen values and vectors if A is square and symmetric"
  286.   (let ((sv (sv-decomp a)))
  287.     (list (coerce (nth 1 sv) 'list) (column-list (nth 0 sv)))))
  288. |#
  289.  
  290. (defun eigenvalues (a) 
  291. "Args: (a)
  292. Returns list of eigenvalues of square, symmetric matrix A"
  293.   (first (eigen a)))
  294.  
  295. (defun eigenvectors (a) 
  296. "Args: (a)
  297. Returns list of eigenvectors of square, symmetric matrix A"
  298.   (second (eigen a)))
  299.  
  300. ;;;
  301. ;;; This is a hack, but it will do for now.
  302. ;;;
  303. (defun eigen (a)
  304. "Args: (a)
  305. Returns list of list of eigenvalues and list of eigenvectors of square,
  306. symmetric matrix A"
  307.   (let* ((n (array-dimension a 0))
  308.          (c (second (chol-decomp a)))
  309.          (sv (sv-decomp (+ a (diagonal (repeat c n))))))
  310.     (list (coerce (- (second sv) c) 'list) (column-list (first sv)))))
  311.  
  312. (defun accumulate (f s)
  313. "Args: (f s)
  314. Accumulates elements of sequence S using binary function F.
  315. (accumulate #'+ x) returns the cumulative sum of x."
  316.   (let* ((result (list (elt s 0)))
  317.          (tail result))
  318.     (flet ((acc (dummy x)
  319.                 (rplacd tail (list (funcall f (first tail) x)))
  320.                 (setf tail (cdr tail))))
  321.       (reduce #'acc s))
  322.     (if (vectorp s) (coerce result 'vector) result)))
  323.  
  324. (defun cumsum (x)
  325. "Args: (x)
  326. Returns the cumulative sum of X."
  327.   (accumulate #'+ x))
  328.  
  329. (defun combine (&rest args) 
  330. "Args (&rest args) 
  331. Returns sequence of elements of all arguments."
  332.   (copy-seq (element-seq args)))
  333.  
  334. (defun lowess (x y &key (f .25) (steps 2) (delta -1) sorted)
  335. "Args: (x y &key (f .25) (steps 2) delta sorted)
  336. Returns (list X YS) with YS the LOWESS fit. F is the fraction of data used for
  337. each point, STEPS is the number of robust iterations. Fits for points within
  338. DELTA of each other are interpolated linearly. If the X values setting SORTED
  339. to T speeds up the computation."
  340.   (let ((x (if sorted x (sort-data x)))
  341.         (y (if sorted y (select y (order x))))
  342.         (delta (if (> delta 0.0) delta (/ (- (max x) (min x)) 50))))
  343.     (list x (|base-lowess| x y f steps delta))))
  344.